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The recent analysis, based on the mean-field approximation (MFA), has predicted that the critical 
quantum collapse of the bosonic wave function, pulled to the center by the inverse-square potential 
in the three-dimensional space, is suppressed by the repulsive cubic nonlinearity in the bosonic 
gas, the collapsing ground state being replaced by a regular one. We demonstrate that a similar 
stabilization acts in a quantum many-body system, beyond the MFA. While the collapse remains 
possible, repulsive two-particle interactions give rise to a metastable gaseous state, which is separated 
by a potential barrier from the collapsing regime. The stability of this state improves with the 
increase of the number of particles. The results are produced by calculations of the variational 
energy, with the help of the Monte Carlo method. 
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I. INTRODUCTION 

It is well known that the attractive inverse-square potential, 

U{r) = -Uor~^, ( 1 ) 

causes the quantum collapse (alias “fall onto the center”) of a wave function of the nonrelativistic particle in three 
dimensions (3D). Potential (1) drives the quantum collapse in its critical form [1]: if Uq exceeds a critical value [which 
is 1/8 in the present notation, see Eq. (2) below], the corresponding Schrddinger equation fails to produce the ground 
state (GS) (or, loosely speaking, in this case the energy of the GS falls to E —>• —oo, while the GS radius collapses 
to zero). The critical form of the quantum collapse is explained by the fact that the quantization breaks the scaling 
invariance of the classical mechanical system with the same potential (1). The effect is also named quantum anomaly, 
or ^^dimensional transmutation” [2]. 

The first solution of the problem of the missing GS was proposed in terms of the linear quantum field theory (rather 
than quantum mechanics, where the problem arises) [2]. That approach actually postulated the existence of the GS, 
with an arbitrary spatial size, which was used as a basic scale for the application of the quantum-field renormalization 
procedure. 

A different solution of the quantum-collapse problem was proposed Refs. [3] and [4], which introduced, in lieu of the 
single quantum particle, a bosonic gas, pulled to the center by potential (1), and described, in the framework of the 
mean-field approximation (MFA), by the Gross-Pitaevskii equation (GPE) for the single-particle wave function [3]: 

^V’ + 47ra^|^pA , (2) 

ot m \ 2 r‘^ J 

where m is the particle’s mass. A key ingredient which secures the GS against the collapse in the framework of the 
GPE is the self-repulsive nonlinear term in Eq. (2). The strength of the self-repulsion is determined by the s-wave 
scattering length, Og, of two-particle collisions [5]. Similarly, the collisional stabilization of a two-component gas pulled 
to the center by potential (1) is predicted by the system of two coupled GPEs [4]. 

The physical realization of the setting, which was proposed in Ref. [3], assumes an ultra-cold gas of molecules 
carrying a permanent electric dipole moment, d (LiGs [6] and KRb [7] molecules may be appropriate for this purpose), 
which are pulled to a point-like electric charge Q placed at the central position (in particular, it may be an ion held 
by an optical trapping potential [8]). In this case, the strength of potential (1) is Uq = \Q\d. An alternative to the 
molecular gas may be provided by a gas of atoms kept in a long-lived Rydberg state [9], which also carry a dipole 
electric moment. The long-range dipole-dipole interactions in the same setting were also taken into account in the 
mean-field theory elaborated in Refs. [3] and [4]. This was done assuming that an individual dipole moment interacts 
with the mean electrostatic field, created by all other particles, i.e., the GPE was augmented by the term accounting 
for the coupling of the dipole moment to the local electrostatic field, and supplemented by the respective Poisson 
equation. As a result, the strength of the contact repulsive nonlinearity was subject to an effective renormalization, 
with an addition induced by the dipole-dipole coefficients: Og —> Ug -I- md^/h?. 
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Thus, it has been demonstrated in Refs. [3] and [4] that the repulsive nonlinearity suppresses the quantum collapse, 
creating the missing GS at all values of Uq, both above and below the above-mentioned critical one, Uq = 1/8. An 
estimate, based on relevant values of the underlying physical parameters, predicted the radius of the restored GS to 
be on the order of a few microns. At Uq > 1/8, a phase transition in the newly built GS was found at Uq = 1/2. At 
that point, the GPE features a breakup of the analyticity in the GS wave function, and a change in the structure of 
its wave function at r —>■ 0 (however, the GS survives the transition). This quantum phase transition resembles those 
found in other systems [12]. 

In addition to Eq. (2) with the isotropic attractive potential (1), its version with the symmetry reduced from 
spherical to cylindrical, which corresponds to potential 17or“^ cos 6, where 6 is the angular coordinate in the 3D space, 
was elaborated in Ref. [10]. Eurthermore, in Ref. [3], the same problem was considered in the two-dimensional (2D) 
geometry with attractive potential (1), for which quantum mechanics gives rise to the collapse for any Uq > 0. Aside 
from the aforementioned realizations, the 2D system models a gas of polarizable atoms that do not carry a permanent 
electric moment, while an effective moment is induced in them by the electric field of a uniformly charged wire piercing 
the system’s plane [11]. Similarly, a magnetic moment may be induced by a current filament set perpendicular to 
the plane [3]. In either case, the 2D version of potential (1) is generated by the attraction of the induced moment to 
the source which induces it. The difference from the 3D setting is that the cubic nonlinearity does not suppress the 
quantum collapse in the 2D setting, only the quintic term (if it is a physically relevant one) being sufficiently strong 
for that [3]. 

The results demonstrating the stabilization of the bosonic gases, against the fall onto the center, by repulsive 
interactions were obtained in Refs. [3], [4], and [10] by means of the MEA. The objective of the present work is to 
address a natural problem suggested by those results, namely, whether similar effects can be found beyond the limits 
of the MFA, i.e., in the framework of a direct analysis of a system of quantum particles pulled to the center by 
potential (1). The corresponding many-body model is formulated in Section II. As a preliminary step, the analysis 
of the single-particle situation is carried out, using the variational approach, in Section III. The mean-field study 
of the many-body system based on the local-density approximation is carried out in Section IV. Physically relevant 
interparticle interaction potentials, namely, of the hard- and soft-core types, are introduced in Section V. The Monte 
Garlo method, used for calculations of the variational energy, is outlined in Section VI. The most essential results, 
obtained by means of this method, are summarized in Section VII, for the many-body system composed of 2 to 
10000 interacting particles. The consideration of a relatively small number of particles, N, is essential, as the MFA 
should be definitely valid for JV ^ 1, in which case the direct analysis of the interactions between particles may only 
produce post-MFA effects, such as two-body correlations. We infer that, while the collapse remains possible in the 
many-body system, the repulsive interactions between particles give rise to a metastable gaseous state, separated 
from the collapsing GS by a sufficiently tall energy barrier. Thus, the possibility of the collisional stabilization of the 
quantum gas attracted to the center by potential (1), predicted by the MFA, is justified by the many-body theory. 
The paper is concluded by Section VIII. In Appendix IX we derive the virial theorem for particles interacting via 
zero-range pseudopotential. 


II. THE HAMILTONIAN 

We study a system of N particles with repulsive interactions, in the presence of the inverse-square potential (I), 
which pulls all the particles to the center. As said above, the system may be unstable against the quantum collapse, 
alias fall onto the center, while the MFA, which takes into regard both the contact and long-range two-particle 
interactions, predicts a possibility of the stabilization [3]. In terms of the quantum-mechanical description, the many- 
body Hamiltonian is 


^^int(|rJ -rfcl), 


( 3 ) 


where Uq is the strength of central potential (1), rj are coordinates of the j-th particle in the 3D space, m is, as in 
Eq. (2), the particle’s mass, and Vint{r) is the repulsive interaction potential between the particles. In the MFA, Vintir) 
is characterized solely by the s-wave scattering length [5], plus the dipole moment which determines the additional 
mean-field interaction of the given particle with the electrostatic field created by all others [3]. Two basic forms of 
the interaction potential chosen for the analysis in the present work are specified below. 

As mentioned above, the attractive potential (1) is a peculiar one, as it scales with distance r in exactly the same 
way as the Laplacian operator of the kinetic energy; as a result, the free-particle problem, with ldnt(|rj — rfc|) = 0, 
has no intrinsic length scale. In other words, the system’s properties scale, following variation of the particle density, 
making phase transitions impossible. This feature is in a stark contrast with the usual situation, where a phase 
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FIG. 1: (Color online) Variational energy (5) for a single particle with 0 = 1 versus (3 for values of the strength of the central 
attractive potential close to the critical value, Uo = 1/8. The upper and lower curves correspond, severally, to free state at 
Uo < 1/8, and a collapsed one at Uo > 1/8. 


transition may be reached by altering the density. The absence of a characteristic length implies that the state of 
the system is fully defined by the attraction strength, [/q, and there is the above-mentioned critical value, Uq = 1/8, 
above which the collapse occurs in the quantum-mechanical system [1]. The terms in Hamiltonian (3) which account 
for the attraction to the center and kinetic energy scale linearly with the number of particles, N, while the energy of 
the repulsive interparticle interactions scales as 7V^. This fact suggests the possibility of the stabilization against the 
collapse in the framework of the many-body theory. 

One may foresee two opposite scenarios related to the possibility of the collapse in the system. In the collapsed 
state, particles fall onto the center, rj —>• 0, creating a diverging density, hence the total energy is negative and diverges 
too. On the other side, in a gaseous state the particles are not bound to the center and fly away, —>■ oo, leading 
to an asymptotically zero density and zero energy. To test what scenario is energetically preferable, we will add an 
artificial Gaussian factor to the wave function, which will localize particles within a finite distance from the center. 

Thus, the single-particle wave function is chosen in the form of a product of a Gaussian and a power-law radial 
term. 


/i(r) = r^exp(-ar2), (4) 

where the inverse Gaussian width, a > 0, determines the inverse localization length. It affects the system’s size and, 
consequently, the density. Alternatively, parameter a can be interpreted in terms of an effective external harmonic 
conhnement with frequency w = 2ahlm. On the other hand, the shape of the cloud at Vj —> 0 is controlled by 
parameter P in ansatz (4). 


III. THE SINGLE-PARTICLE SOLUTION 


The single-particle problem, dehned by Hamiltonian (3) with N = 1, can be studied by means of the variational 
method, treating a and /3 in ansatz (4) as variational parameters. In the single-particle sector, the system is steered 
by the competition of the external potential and kinetic energy, while the interparticle potential, Vint(|ri — rj|), does 
not appear. The variational energy, Mar = [/ /i(r) dr] ^ f /i(r)id/i(r)dr, with fi taken as per Eq. (4), can be 
evaluated explicitly: 


E^^'> = a 


8Uo-l 
2(1-f 2/3) 


( 5 ) 


For a fixed localization size, a = const, this energy is a decreasing function of /3 if Uq is smaller than the critical value, 
Uo = 1/8 (the same which was mentioned above), and an increasing function for Uq > 1/8, see Fig. 1. By varying /3, 
we probe how the energy responds to the increase or decrease of the density at the center of the cloud. 

For Uq < 1/8, the minimum of energy (5) with respect to /3 is A = a, which is reached at /3 —>• oo, corresponding 
to a delocalized state with the particles flying away from the center. The minimum of the energy with respect to a 
is F = 0, which is attained at a = 0, i.e., for zero density, when the localization factor in Eq. (4) is absent. 
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The denominator in expression (5) vanishes at /3 = —1/2, when the energy diverges. Close to this point, at 
P = —1/2 + e with small e, the energy is (1/8 — Uq + e/2). Depending of the sign of 1/8 — C/q, the 

energy diverges to negative or positive infinity. Thus, for Dq > 1/8, the collapsed state is realized, corresponding 
to a —^ oo (zero localization radius). On the other hand, for /3 —^ —1/2 in the absence of the confinement, a = 0, 
the density profile at small distances is p{r) oc //(r) = = 1/r, the respective expression for the potential energy, 

/p p{r)Uo/r^ Airr^dr, featuring a logarithmic divergence. 

The divergence at Uq > 1/8 will persist in the many-body system. This can be seen as follows: one can place 
all particles but one off the center, which will cost a finite energy. Still, at Uq > 1/8, the part of the total energy 
corresponding to the single particle staying near the center will be divergent and negative. That single particle falls 
onto the center. One realizes that if the potential energy of the pair-interaction remains finite, all other particles 
will also fall onto the center. In other words, the absolute stability threshold remains the same, Uq = 1/8, also in 
a many-body system. On the other hand, a metastable state may appear in the many-body system with repulsive 
interparticle interactions, as predicted by the mean-field consideration [3]. 


IV. THE MEAN-FIELD THEORY BASED ON THE LOCAL-DENSITY APPROXIMATION 


Here we address the problem of stability of the many-body system within the framework of the local-density 
approximation (LDA). This method is commonly applied to many-body systems in the presence of an external 
potential, when the chemical potential of the homogeneous system is known. We use the mean-field expression for 
the chemical potential, /ihom = gn , where g = Anh^as/m is the coupling constant. This choice corresponds to the 
short-range interaction potential determined by the s -wave scattering length as per the lowest Born approximation. 
Mathematically, pseudopotential V{r) = g5{r) leads to the same result. It is important to note that this is a soft-core 
interaction, as it allows particles to overlap. The pseudopotential interaction permits to derive the virial theorem in 
the presence of the central potential, see Eq. (26) in Appendix. It is of interest to introduce inhomogeneity by an 
external trap and check if the energy remains positive or negative once the trap is suddenly removed. 

In the framework of the LDA, the chemical potential in the presence of an external potential is approximated by 
the sum of the local chemical potential pioc = gn, where, this time, n is a function of the coordinates, rather than a 
constant, and the external potential, 


1^0 1 2 2 

P = Pioc -y + ttWo; r , 

2 


( 6 ) 


where the harmonic-oscillator confinement, with the respective length scale, Oho = ybJ(rnLo), is added to make the 
size of the system finite. Solving Eq. (6) for the density, one obtains the following density profile: 


1 j p — + Uqt at r < Atf, 

1 0, at > i?TF, 


( 7 ) 


where the cloud’s radius is given by the Thomas-Fermi (TF) expression, i?TF = \JP+ \/P^ E 2mUouj'^/{^JrnuJ), which 
diverges for a shallow trap, at w —>■ 0. Close to r = i?TF, the density profile follows a typical inverted-parabola shape 
of a trapped Bose condensate. The density at the center features an integrable divergency reflecting the presence of 
the attractive central potential. The value of the chemical potential itself is fixed by the normalization condition, 
47r n(r)r^dr = N. 

The system’s energy is obtained by integrating the chemical potential: E = p{N')dN' . We are interested to 
find the release energy, which remains in the system when the harmonic trapping is suddenly switched off, Erei = 
E — \N mu)'^{r'^), where (...) stands for the average value. It is especially important to find out if is positive or 
negative. 

In the limit of the vanishing central potential, Uq = 0, the size of the cloud is defined by the balance between the 
chemical potential (which depends on the particle-particle interaction strength) and the harmonic confinement [5], 

i?TF = \/2p/muj^. The respective chemical potential, p = \hjjj (ISiVos/aho)^’^^, is always positive, hence the respec¬ 
tive release energy, Arei = {2/7)pN, is positive too. 

In the opposite limit of strong attraction to the center, the chemical potential, p ~ /(9V^/i^a^), is negative 

and independent of the trapping frequency, while the TF radius is inversely proportional to the strength of the 
attractive potential: Atf = iNh^Os/ {2mUo). In this case, the system’s energy diverges for on the lower integration 
bound corresponding to small system’s sizes, E —>■ —oo. Thus we conclude, that, in the framework of the LDA, the 
collapse is always possible, as far as the energy is considered. 
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It is important to note that the kinetic energy, which is a key ingredient precluding the onset of the collapse, is not 
fully taken into account by the LDA. A more complete mean-field Gross-Pitaevskii theory, which takes into account 
the kinetic energy properly, predicts the suppression of the collapse [3]. Thus, it is obviously interesting to see how 
beyond-mean-field effects modify that result. 


V. SOFT- AND HARD-CORE INTERPARTICLE POTENTIALS 


In the single-particle problem, the collapsed state results in a density profile which is a delta-function at the center, 
p(r) = (5(r). Then, two essential questions arise: (i) whether a similar singular state, p(r) = N6{r), corresponds to 
the GS of a many-body system containing N particles; and (ii) if a nonsingular gaseous state may be metastable, i.e., 
dynamically stable against small perturbations, even if the GS corresponds to the collapse. 

The answer to question (i) depends on the form of the interparticle potential, Vint(r’), at small r. If the interaction 
features a soft core, i.e., there is an upper bound on the interparticle interaction potential Vq, such that V{r) < Vq at 
r —>■ 0, then there also exists an upper bound on the particle-particle interaction energy: E-mt < VoN{N — l)/2. At 
the same time, the potential energy corresponding to the central attractive potential diverges to —oo at Uq > 1/8. 
As a result, any soft-core potential is not sufficient to prevent the system from the transition to the collapsed GS (fall 
onto the center). On the other hand, for physically realistic interparticle potentials with a hard core (for example, 
of the Van der Waals type), with Vg oo, the fully-collapsed state, with p{r) = N6{r), would lead to a diverging 
positive contribution to the energy, which may prevent this state from being the GS. 

The answer to question (ii) depends on the number of particles, N, and parameters of the interaction. It may 
happen that, while the GS corresponds to the fall onto the center, there is a barrier between the collapsed and 
nonsingular gaseous states, making the latter one metastable. Indeed, the interparticle-interaction energy, Eint, 
depends quadratically on N, while the energy of the attraction to the central potential scales linearly with N, which 
suggests that the systems may effectively stabilize itself for N large enough. 

To study the expected scenarios for the soft- and hard-core potentials, we consider two respective models for Vint(?'): 
the hard-sphere potential of diameter R, 


and its soft-sphere counterpart. 


Thard {r ) 


OO, r < R 
0, r>R’ 


Vsoft (r) 


Vo, r < R 

0, r > i? ’ 


( 8 ) 

(9) 


with finite Vq in the latter case. By varying height Vq of the soft-sphere potential one can change the value of the 
respective s-wave scattering length, which is 


Os = R[1 — tanh{kR)/{kR)], (10) 

where k = y/mVojh is the momentum associated with the height of the soft sphere. For hard-sphere potential (8), the 
effective s-wave scattering length is identical to the diameter of the sphere, as = R. Note that the s-wave scattering 
length (8), corresponding to the soft-sphere potential, is smaller than the diameter, approaching it in the limit of 
Vb oo. As mentioned above, for dilute gases it is expected that Ug is the single parameter which determines effects 
of the interparticle interactions [5]. Below, we compare the interaction potentials with the same s-wave scattering 
length Us but different diameters. 


VI. THE MONTE CARLO METHOD 

An efficient way to calculate the energy of a many-body system is to use the Monte Carlo technique. We resort to 
the variational Monte Carlo method, which samples the probability distribution, p = for a known many-body 
wave function, tp, allowing one to obtain the variational energy as a function of trial parameters, such as a and /3 in 
Eq. (4). The standard Metropolis algorithm [13] is used for the implementation of the method. 

The many-body trial wave function is chosen as a product of single-particle terms, fi(r), and a pairwise product 
of two-particle Jastrow terms [14], f 2 {r): 

N N 

iIi(yi,...,vn) = n-^i(’^*)n-^2(|ri-rjl) 

i=i i<j 


( 11 ) 
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The single-particle factor is taken in the same form as the above single-particle ansatz, given by Eq. (4), with parameter 
a controlling the degree of the particle’s localization. In particular, a —>■ 0 corresponds to the limit of free delocalized 
particles, while a —>■ oo implies the fall onto the center. The Jastrow factor f 2 {r) in Eq. (11) is chosen as a solution 
of the two-body scattering Schrbdinger equation. In this way, the short-range physics in taken into account correctly, 
and interparticle correlations, which are important in the context of the metastability of a many-body system, are 
introduced. For the hard-sphere potential, the two-body solution is given by 




0, r < R 

1 — R/r^ r > R 


( 12 ) 


Note that, for distances larger than the hard-sphere’s radius, solution 1 — i?/r is the same as for a pseudopotential 
with the respective s-wave scattering length, Ug = R. 

For the soft-sphere potential (9), the two-body scattering solution is [1] 




Asinh(fcr)/r, r < R 
1 — ttg/r, r > R 


(13) 


where the momentum scale is fc = ^JmVo/h, see Eq. (10), and constant A is determined by the condition of the 
continuity of / 2 (t’) at r = R. Two hard spheres cannot overlap, hence the Jastrow factor (13) vanishes at r < i?, as 
written in Eq. (12). On the other hand, for the soft spheres, there is a finite possibility of finding two particles at the 
same place, which allows the Jastrow factor to remain finite at r = 0: = 0) = Ak r; (i? — as)kexp{—kR) > 0, 

with Qs < R. The possibility for two particles to overlap gets exponentially suppressed as the height of the soft sphere 
increases. Two-body scattering solutions (12) and (13) coincide in the limit of Vb — oo. 


VII. NUMERICAL RESULTS FOR THE MANY-BODY SYSTEM 
A. The two-particle system 

First, we address the response of the two-particle system to compression. To this end, we calculate, by means of 
the Monte Carlo method, the dependence of the variational energy on the density, which is controlled by localization 
parameter a in Eq. (4). Figure 2 shows the variational energy for C/q = 1 and /? = 0. It is seen that pairs of particles 
interacting via the hard-sphere potential are always stable against the many-body collapse. In the system with zero 
density, a = 0, the energy is also zero, E = 0. There is a minimum in the energy at a finite value of a, meaning 
that the GS corresponds to a bound state {E < 0) of a finite size. For large values of a (i.e., for a small size of the 
localized system), the energy becomes large due to the strong interparticle repulsion. A soft-sphere potential with a 
small size, R = l.los, which, according to Eq. (10), means large potential height Vq, follows behavior similar to that 
of the hard-sphere potential for the values of a shown in Fig. 2. 

For larger sizes of the soft sphere and smaller potential strengths, Vb, a second minimum in the variational energy 
becomes visible at finite values of a reported in Fig. 2. The absolute energy minimum stays at a —> oo, implying, as 
said above, that the GS corresponds to the fall onto the center (quantum collapse). The gaseous state is, nevertheless, 
metastable if there is an energy barrier separating the two minima. In the present case, we observe that a finite 
barrier is present at Uq < 1.4. For a stronger attraction to the center, there is only one minimum, corresponding to 
the collapse. 

Next, we aim to approximate the energy of a highly compressed state. For a very strong localization, a —>■ oo, the 
bosonic particles stay on top of each other, hence the potential energy of the two-body interaction is well approximated 
by the height of the soft-sphere (9), Ei^t = Vq- The two-body repulsion leads to an additional increase of the system’s 
size, thus decreasing the kinetic energy per particle, ifkin/-N, and reducing the attraction to the center, Ecentred/N. 
We find numerically that the net effect amounts to an increase the sum of both energies by Vb, against the one- 
particle value (5), Fkm + ^center = 2E^^'> + Vq. As a result, the total energy of the two-particle system, E^'^\ can be 
approximated by 


£:(2) ^ -f 2Vb. (14) 

The corresponding asymptotic behavior is shown in Fig. 2 by dashed lines. In the limit of a —>■ oo, the energy of the 
two-particle system with the soft-sphere interaction is divergent and negative, corresponding to the quantum collapse. 




FIG. 2: (Color online) The variational energy per particle as a function of the inverse-Gaussian-width parameter, a, for Uq = 1, 
Os = 0.1 and N = 2 particles, with the hard- and soft-sphere potentials. Solid lines pertain to the soft-sphere potential with 
the diameter of the sphere, R, ranging from l.las to 2as (from the upper to the lower curve). The dashed line corresponds to 
the hard-sphere potential (i? = Os). The dash-dotted line depicts the single-particle energy given by Eq. (5). 



FIG. 3: (Golor online) The energy per particle for the soft-sphere interaction potential as a function of the inverse-Gaussian- 
width parameter, a, for Uo — 1, aa = 0.1, R = 1.3as and the number of particles N — 2; 3; 4; 5; 10,100,1000, 10000 (larger 
number of particles has larger value at the maximum). Solid lines, variational result; dashed lines, the asymptotic energy of 
the fully-collapsed state, as per Eq. (16); dash-dotted line, typical energy associated with the Gaussian localization, given by 
Eq. (15). 


B. Many-body systems 

It has been demonstrated above that the energy barrier may be present between the gaseous and collapsed states 
for the two-particle system, N = 2. With the increase of N, the interparticle interactions improve the stability of the 
former state. 

Figure 3 shows the variation energy, calculated by the Monte Carlo method for a fixed radius of the soft sphere, 
R = l.Sos, for a wide range of the number of particles, N = 2 — 10000. For localization parameter a = 0 the density 
vanishes and the energy is equal to zero. For small values of a, which corresponds to the weak localization, the energy 
might be negative (not visible on the log-log plot of Fig. 3). As the localization gets tighter, the energy becomes 
positive, as the two-body interaction helps the system to resist the compression. For very tight localization, a —>■ oo, 
the fall onto the center is observed with the energy diverging to minus inhnity. 

For large N, the energy calculated with ansatz (4) does not immediately lead to the fully collapsed state. The 
localization energy, proportional to hui of the harmonic trapping, may become a dominating term in the energy, while 
the system’s size is still large enough, so that the fully-collapsed state is not realized. The energy in the corresponding 
regime is numerically approximated as 


E = NCa 


(15) 


with (7 = 4, and is shown in Fig. 3 by the dashed-dotted line. For still tighter localization, it has been found that 
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FIG. 4: The energy barrier between the state with a = 0 and a —>■ oo as a function of the number of particles, N, for the data 
shown in Fig. 3. The dashed line depicts the asymptotic approximation (17) for the large system. 


the following expression well approximates the energy of the system in the limit of the full collapse, when all particles 
overlap: 

E = + N{N-l)Vo. (16) 

Similarly to the case of the two-body system [see Eq. (14)], the net effect of the interparticle interactions, revealed by 
our calculations, is two times larger than the potential energy, Ei^t = VoN{N — l)/2. The asymptotic energy (16) is 
shown in Fig. 3 by dashed lines. 

A clear conclusion is that the increase of the number of particles indeed causes strong rise of the potential barrier 
which stabilizes the energy minimum corresponding to the gaseous state. This can be also concluded from Eq. (16), 
where the contribution due to the repulsive interactions scales as N'^ for large N, while the term corresponding the 
attractive central potential scales as N. 

The energy barrier between the state described by ansatz (4) and the free state with zero energy, Ebarrier, is 
estimated as the maximum value of the energy per particle (see Fig. 3), and is shown in Fig. 4. For a large system 
size, the barrier can be approximated by comparing the two basic energy scales given by Eqs. (15) and (16). This 
resulting asymptotic approximation is 


Abarrier = 6NVo/{8Uo - 3 + 2C), 
which is shown in Fig. 4 by the dashed line. 


(17) 


VIII. CONCLUSIONS 

This work is stimulated by the recent results showing that, in the framework of the MFA (mean-field approximation), 
the quantum collapse (alias “fall onto the center”) of bosons pulled by 3D attractive potential (1) can be suppressed 
by the repulsive collisional nonlinearity in the ultracold bosonic gas, which is described by the GPE (Gross-Pitaevskii 
equation) [3, 4]. Accordingly, the collapsing GS (ground state) is replaced by a regular one with a finite density. Here, 
going beyond the framework of the MFA, we demonstrate that this prediction finds its counterpart in the many-body 
system governed by the same Hamiltonian: while the collapse is still possible in the system, the repulsive interactions 
create a non-singular metastable gaseous state, which is isolated from the collapsing one by a tall potential barrier. 
The metastability quickly enhances with the increase of the number of bosons, N, starting from V = 2, in agreement 
with the expectation that the MFA should be valid for N large enough. In the same time, the absolute stability 
threshold remains the same, in terms of Eq. (1), Uq = 1/8, as for the single particle [Ij. The results have been 
obtained through calculations of the variational energy, implemented by means of the Monte Garlo technique. 

The analysis reported here may be naturally expanded in other directions. In particular, it may be possible to 
construct stabilized states carrying the angular momentum, i.e., counterparts of vortical modes, predicted in Ref. [3] 
with the help of a crude version of the MFA. Further, following the MFA results for the 2D geometry reported in 
Ref. [3], it may be interesting to consider the many-body system in the 2D setting. The conclusion that only the 
quintic term may suppress the quantum collapse in that case, suggests that the stabilization of the 2D state may 
demand adding three-body interactions to the many-particle Hamiltonian. Lastly, a challenging problem is to explore 
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a similar possibility for fermions pulled to the center by potential (1). The simplest approach may then be to consider 
a balanced Fermi gas, with repulsive interactions between equal numbers of spin-up and spin-down particles [15, 16]. 

B.A.M. appreciates hospitality of ICFO - Institut de Ciencies Fotoniques (Barcelona, Spain). G.E.A. acknowledges 
partial financial support from the MICINN (Spain) Grant No. FIS2014-56257-C2-1-P. The Barcelona Supercomputing 
Genter (The Spanish National Supercomputing Center - Centro Nacional de Supercomputacion) is acknowledged for 
the providing access to computational facilities. 


IX. APPENDIX: THE VIRIAL THEOREM 


In this Appendix we derive the virial theorem for the pseudopotential interaction in presence of both the central at¬ 
tractive potential (1) and harmonic-oscillator trapping. For the sake of completeness we consider arbitrary dimension, 
D = 1,2,3. 

In the general form, the Hamiltonian is given the following sum 


where 


77 — Ekin “f E/central “t” ^ho T ^int: 


^ 2 


(18) 


(19) 


i=i 


is the kinetic energy, 


N 


Uo 


^central — '^ 2 


describes the attraction to the central potential, 


r • 

3 = 1 


N 


( 20 ) 


Eho = ^ 


( 21 ) 


1=1 

is the potential energy of a the harmonic-oscillator confinement, and 

N 

Eint = ~ I’feL 

j<k 


( 22 ) 


determines the interparticle interactions mediated by the pseudopotential. Strictly speaking, in Ajnt regularization 
of the ^-functional interaction is needed in three dimensions, but it will not affect the scaling argument used here, 
therefore we disregard the regularization. 

We now aim to check how the rescaling of particles’ coordinates, r —>■ ^r, affects averages of the different terms in 
the Hamiltonian. The normalized GS wave function transforms according to 




...,^rjv) 




\// dri...drN 

while the energy becomes 

E= f dr,...drN = e{Ekin) + e"(4ental) + r"(^ho) + 


(23) 


(24) 


The principle that the energy must be minimal for the GS imposes the following condition for a small variation 
around the ground-state solution with = 1: 


JE 


= 0 . 


«=i 


(25) 
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As a result, we obtain the following relation (recall that d is the spatial dimension): 

2(i?kin) + 2(i?central) ~ 2(^ho) + D {Eint) = 0, (26) 

which is the virial theorem for Hamiltonian (18). In the absence of central potential (1), Eq. (26) amounts to the 
known result [for D = 3, see Eq. (44) in Ref. [17]]. 
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